Black-hole threshold solutions in stiff fluid collapse 
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Numerical studies of the gravitational collapse of a stiff {P — p) fluid have found the now familiar 
critical phenomena, namely scaling of the black hole mass with a critical exponent and continuous 
self-similarity at the threshold of black hole formation. Using the equivalence of an irrotational stiff 
fluid to a massless scalar fleld, we construct the critical solution as a scalar field solution by making 
a self-similarity ansatz. We find evidence that this solution has exactly one growing perturbation 
mode; both the mode and the critical exponent, 7 ~ 0.94, derived from its eigenvalue agree with 
^ , those measured in perfect fluid collapse simulations. We explain why this solution is seen as a critical 

\^ • solution in stiff fluid collapse but not in scalar fleld collapse, and conversely why the scalar field 

' critical solution is not seen in stiff fluid collapse, even though the two systems are locally equivalent. 

o ■ 
l> 
o 

o 

■ I. INTRODUCTION 

I I Solutions of Einstein's equations that lie at the threshold of black-hole formation have provided a unique insight 
bJO into gravitational dynamics during the past decade. The advent of sophisticated numerical techniques and powerful 
^ ^ computers has facilitated the study of this regime. 
• rH , Critical phenomena were first discovered in the gravitational collapse of a real, massless scalar field, in the evolution 
^\ of generic 1-parameter families of asymptotically flat initial data During the time evolution, the matter either 
disperses (subcritical data) or collapses to form a black hole (supercritical data). The critical solution exists precisely 
at the threshold of black hole formation. This solution has several interesting properties which also are at least 
approximately shared by near-critical solutions. First, the masses of black holes formed in supercritical evolutions 
obey a universal scaling law, 

MBH«b-p*P, (1.1) 

where p is the parameter of a given 1-parameter family of data, and p* is its critical value such that black holes 
are formed for p > p*. Thus, in this model, infinitesimally small black holes can be formed by adjusting p to be 
sufficiently close to p* . Second, near-critical evolutions go through a universal intermediate phase that exhibits self- 
similar echoing near the origin. The echoes are periodic in logarithmic time, so that the critical solution (universal 
intermediate attractor) is, in fact, discretely self-similar. Empirically, it is found that the critical exponent, 7 ~ 0.374, 
is the same for all 1-parameter families of massless scalar field initial data. 
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As they have now been found in a wide variety of general relativistic self gravitating systems, critical phenomena 
are considered to be generic features of "tuned" gravitational collapse. Interested readers are directed to for 
a comprehensive review. In brief, the threshold of black hole formation singles out well-defined critical solutions 
which typically have additional symmetry (beyond that which might have been imposed via the formulation of the 
model under consideration), and which, though unstable by construction, tend to have precisely one unstable mode in 
perturbation theory. The critical solutions are typically unique (up to a dynamically irrelevant overall rescaling): the 
unstable mode is consequently also unique, as is the growth factor associated with the mode. The growth factor can 
then be immediately related, via ideas familiar from renormalization group studies, to a scaling exponent defined, for 
example, via (|1T|) |||,||. 



The critical solutions which have been observed so far have either a time-translational or a scale-translational 
symmetry [^,^; these two types have also been dubbed Type I and Type II, respectively. Within each of the 
two symmetry categories, solutions can be further differentiated according to whether the extra symmetry appears 
continuously or discretely. Thus, Type I solutions are static or periodic, and display a mass gap at threshold. As 
one tunes closer and closer to criticality, the static/periodic intermediate attractor (star-like solution) persists for a 
lifetime r ~ cr In |p — p*| where a is the reciprocal Lyapunov exponent associated with the unstable mode. Type II 
solutions, on the other hand, are either continuously or discretely self-similar (CSS or DSS), exhibit infinitesimal mass 
at threshold, and generically have naked singularities at t he o rigin in the precisely critical limit. In the supercritical 



case, the black hole mass scales according to the relation (1.1), where again, the scaling exponent 7 is the reciprocal 



Lyapunov exponent of the critical solution's single unstable mode. The remainder of this paper concerns only Type 
II behavior, and the key distinction between the two solutions which will be discussed is that one is CSS while the 
other is DSS. We also note that a key feature of the Type II solutions which have been computed thus far via direct 
tuning of PDE solutions is that, virtually by definition of existing at the threshold of black hole formation, they do 
not contain apparent horizons. 

One of the models in which critical phenomena have been studied most extensively is that of spherically symmetric 
collapse of perfect fluids with the one-parameter scale-free equation of state (EOS) 

P^kp. (1.2) 

Here P is the fiuid's isotropic pressure, p the energy density, and < fc < 1 is a constant (the single parameter of 
the EOS). Critical collapse in this context was first considered by Evans and Coleman for the specific case k = 1/3 
(radiation equation of state) These authors not only found a continuously self-similar (CSS) critical solution via 
evolution of families of initial data, but also constructed that same critical solution more directly by adopting a CSS 
ansatz, and then solving the resulting boundary-value eigenproblem. This work was quickly followed by examination 
of the stability properties of this and related CSS fiuid solutions within the context of perturbation theory . 
Working from the CSS ansatz, these studies identified single-mode-unstablc CSS solutions for < fc < 0.89, indicating 
that the black hole threshold solutions were CSS (and thus Type II) for fc in this range. The absence of solutions 
for k > 0.89 was intriguing, and hinted that the critical solutions for such cases might have different properties — in 
particular, there were suspicions that one or both of the assumptions of 1) regularity at the origin, and 2) continuous 
self-similarity might break down at fc ~ 0.89. However, CSS critical solutions with regular origins for 0.89 < fc < 1 
were subsequently found, through fiuid evolutions PJTo[|, as well as via adoption of the CSS ansatz with, from the 
point of view of the PDE evolutions, no obviously special changes in solution properties at fc ~ 0.89 ||]. 

The existence of a CSS critical solution in the limiting case of the stiff fluid (fc — 1) was particularly surprising 
to many, as the spherically symmetric stiff fluid is formally equivalent to a spherically symmetric real scalar fleld. 
[The technicalities of this equivalence are detailed below in Section II A for the stiff fluid, and in Appendix |^ for an 



arbitrary barotropic equation of state P = P{p)-] Thus, one might have naively expected that the critical solution 
appearing in one system would automatically arise as the critical solution in the other. However, the DSS solution of 
scalar field collapse and the CSS solution of stiff fluid collapse are distinctly different solutions. 

Two fairly obvious questions then arise: 1) Why is the CSS critical solution found in stiff fluid collapse not observed 
in massless scalar critical field collapse? 2) Why is the DSS scalar fleld critical solution not seen in stiff fluid collapse? 
The answer to the second question is relatively trivial, and has been well known in the critical-collapse community 
for some time — the DSS scalar fleld critical solution corresponds to a perfect fluid with negative (as well as positive) 
pressures and densities, and hence cannot arise in fluid collapse scenarios where energies and densities, are, by flat, 
constrained to be positive. As discussed in more detail below, the resolution of the second question is somewhat 
more subtle, but in essence amounts to the observation that the CSS critical solution, when interpreted in the scalar 
fleld context, contains an apparent horizon and therefore cannot and does not sit at the black hole threshold. This 
appears closely related to the phenomena observed in Lechner et al's study of general relativistic non-linear sigma 
models where the absence of CSS solutions at black hole threshold above a certain strength in coupling parameter 
is correlated with the fact that the CSS solutions at those coupling values have apparent horizons. 
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II. A CONTINUOUSLY SELF- SIMILAR (CSS) SCALAR FIELD SOLUTION WITH ONE GROWING 

MODE 



In this section we construct the critical sohition for spherically symmetric stiff fluid collapse using the equivalence 
between the scalar field and the stiff fluid to work in the scalar field variables. We first describe this equivalence 
between the two matter models, give the field equations in Bondi coordinates, then restrict them to a continuously 
self-similar ansatz. We then study the spacetime structure and perturbation modes for two solutions, one of which is 
the desired critical solution. 



A. The correspondence between the stiff fluid and the scalar field 

A stiff fiuid is a perfect fluid with equation of state, P — p. The stress-energy tensor for the stiff fluid can therefore 
be written as 

T"^ ^ p{2u''v!' + g"'') . (2.1) 

Here p is the total energy density of the fluid as measured by an observer moving with a fluid clement, and hence 
having a normalized 4- velocity u". If the 4- velocity is irrotational (V = 0), we can write in terms of the 
gradient of a scalar field. (In spherical symmetry the fluid is automatically irrotational.) We introduce a scalar field 
-0 with future-pointing timelike gradient such that 

Va^=y/2'pua, (2.2) 

where Vq is the metric-compatible covariant derivative. This defines ip up to an additive constant, whose value is 
irrelevant. The inverse relation between the fluid and the scalar field is then given by 

p - -iVcV^V^i/;, (2.3) 



(2.4) 



Expressing the stress-energy tensor (2.1) in terms of tp gives 

yafc ^ v'^VVV - i(VcV'V'=V)ff°^- (2.5) 

The conservation of energy-momentum, WaT'^^ — 0, implies that ifj satisfies the massless, minimally coupled scalar 
wave equation 

VcV^V = 0. (2.6) 

This establishes a local, one-to-one relationship between irrotational (in particular, spherically symmetric) stiff fluid 
solutions and massless scalar field solutions with timelike gradient. The existence of massless scalar field solutions 
that also have spacelike and null gradients is relevant later when we discuss the relation between interpolating families 
of stiff fiuid solutions and those of scalar field collapse. A more general treatment, applicable for equations of state of 
the form P ~ P{p), is given in Appendix 

Although the perfect fiuid and scalar field models are related by this mathematical equivalence, differences remain in 
their physical interpretation. Generally, a fluid is a phenomenological model derived via thermodynamic considerations 
of a large number of particles interacting through elastic collisions. The fluid density, p, and four- velocity, are 
continuum variables deflned in terms of the fundamental discrete variables in a limiting procedure. The fluid energy 
density, for example, is defined in terms of the total energy in an arbitrary volume, in the limit that the volume element 
becomes infinitesimally small. Thus, a fluid is a continuum model of a discrete system containing a thermodynamically 
significant number of particles; the model is expected to fail as the large-particle-number limit is violated, p 0. 
Indeed, the fluid equations become singular in this limit. 

Given the thermodynamic motivation for perfect fluid models, one conventionally imposes physical constraints on 
the fluid's energy density and four-velocity, namely that p > 0, and that m° is timelike. The fiuid equations may 
admit solutions in violation of these constraints, in which case the solution may be discarded, or matched with other 
solutions to create a physical spacetime. One could argue against accepting these physical constraints, but one would 
then be adopting a view of a fluid seemingly at odds with the usual assumption of an underlying discrete system. 

The scalar field, on the other hand, is fundamentally a continuum model, and does not directly represent the 
behavior of a discrete system of particles. Consequently, its solutions are generally not subject to the constraints 
discussed in the previous paragraph. Thus, the assumptions made in constructing the perfect fluid model allow only 
a restricted class of scalar field solutions — those with V^ipVatp < — to be interpreted as physical fluid solutions. 
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B. Coordinate system and field equations 



The general spherically symmetric line clement can be written in Bondi coordinates as 

ds^ = -ggdu^ - 2gdudr + r^{de^ + sin^ edc/)^) , (2.7) 

where g = g{u, r) and 'g = g{u, r). The coordinate r has geometrical meaning in terms of the proper area of r = const. 
2-spheres. We demand that the critical solution be regular at the origin; elementary flatness at r = then dictates 
that g{u, 0) — g{u, 0). The relation between the fluid density and scalar field in these coordinates is 

p^^{2dui^-gdr^) . (2.8) 

A local mass function m{u, r) is defined by 

l-'-^^^g^^V.rV.r^l, (2.9) 
r g 

where r is understood to be a function on spacetime. 

Looking ahead to the application of the CSS ansatz, it is convenient to express the Einstein equations using variables 
adapted to a scale symmetry: 

X = -r/u , T = - In(-u) . (2.10) 

This coordinate change is applied directly to the Einstein equations (not to the line element), which then become 

xg'^g-g, (2.11) 
xilng)' = Anf , (2.12) 

(g/g)- + x{g/g)' = + V-) [{x - g)j + x^] /<? . (2.13) 

Here a prime (') represents differentiation with respect to x, a dot (') means differentiation with respect to r, and we 
have also introduced an auxiliary field, j(u,r), defined by 

x^'=i. (2.14) 



Finally, the scalar wave equation of motion (2.6) becomes 



(V' + ]y + {x- g/2)f = A(g _ 2x) . (2.15) 

Now that the Einstein equations have been written in terms of the new variables, the CSS ansatz will reduce them to 
a system of ODEs. 

C. Continuously self-similar solutions 

Self-similar solutions representing scalar field collapse have been studied by Brady ||ll| and Christodoulou jl^ . We 
present here only those details that are essential to our discussion; the interested reader can consult the above-cited 
references for further details. 

The existence of a homothetic symmetry (continuous scale symmetry) in a spherical spacetime implies that the 
metric coefficients g and g depend only on a; = —r/u, and that the scalar field is of the form 

^p = h{x)+KT, (2.16) 



where h is a function to be determined, and k is a positive constant. Restricting equations (2.11)-(2.15) to the 
self-similar ansatz gives 
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FIG. 1. A spacetime diagram showing the self-similar scalar field solution with ^-ki^ ~ 0.577. The apparent horizon, the 
similarity horizon at x = 1, and the surface on which the energy density p vanishes are all shown. The shaded region, where 
the scalar field gradient is timelike, can be identified with a physical stiff fluid solution. 



xg'^9~l, (2.17) 
x(ln,g)'=4^j2, (2.18) 
g - 5 - ^-kXIk^x - {g-2x){f + 2kj)] , (2.19) 

xTi = j , (2.20) 
x{2x - g)j' = -2kx + j{g - 2x) . (2.21) 

These equations are singular when g = 2x. This corresponds to a similarity horizon in spacetime, that is a nuU 
hypersurface of constant x coordinate. With a slight abuse of terminology, it is also a sonic surface of the corresponding 
stiff fluid. We rescale the coordinate w by a constant factor such that x = 1 when 5 — 2x. Note that this means that 
the usual convention g = g = 1 at the origin x = does not hold. 

In order to be regular in the past, say at u = —1, the threshold solution must be analytic at the origin and at 
the similarity horizon x — 1. This is impossible for general values of k, but we expect to find analytic solutions for 
isolated values of kJ 3| . 

Christodoulou [|l2rhas derived a closed form solution, satisfying these conditions, when Attk^ = 1/3. This is a 
Friedmann-Robertson- Walker cosmological solution, and is not the similarity solution observed in numerical simu- 
lations of gravitational collapse. We use it here only to test our numerical methods on an exact solution. Using 
a two point shooting method to search for analytic solutions at other values of k, we find a second solution when 



Attk^ ~ 0.577. This solution does correspond to the stiff fluid critical solution, as discussed below in Section IIE, but 
first we discuss its properties as a massless scalar field solution. 

The spacetime structure of the Attk'^ ~ 0.577 solution is sketched in Fig. ^ There is a strong curvature singularity 
at r = and u = 0; this is a point, in the sense that it can be surrounded by an arbitrarily small sphere in spacetime. 
The center of spherical symmetry, x = 0, for u < is made regular by our ansatz, and so is the past light cone of 
the singularity, x = 1. The scalar field gradient is timelike for < x < xi, with xi ~ 1.07. In that region, the scalar 
field matter can also be interpreted as a stiff fluid. The surface x — xi is a regular spacelike surface to the future 
of X = 1, on which the scalar field gradient is null. In the fiuid interpretation this corresponds to the limit in which 
the fluid 4- velocity becomes null, and the comoving fluid energy density goes to zero. For x > xi, the scalar field 
gradient is spacelike, and the scalar field cannot be interpreted as a stiff fluid given the usual physical constraints on 
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TABLE I. Numerical results for the modes of the continuously self similar scalar field solution discussed in the text. Listed 
in columns 3 through 7 are the mode eigenvalues, A, and, where relevant, corresponding scaling exponents, 7 = 1/A, computed 
using the two methods described in text. Also listed is the scaling exponent 7pde estimated from solution of the full equations 
of motion Quoted uncertainties are estimates based on the maximum discontinuity at the fitting point for the shooting 
method and convergence tests for the matrix method. 



Mode 


Acxact 


'^shoot 


'Yshoot 


'^matrix 


'Ymatrix 


7PDE 


physical unstable 




1.0654 ± 0.0005 


0.9386 ± 0.0005 


1.055 ± 0.01 


0.95 ± 0.01 


< 0.96 


gauge shift 


1 


1.0000 ± 0.0005 




1.013 ± 0.01 






gauge rescale 





0.0000 ± 0.0005 




0.000 ± 0.00 







the scalar field's gradient. The spacehke surface x = X2 — 1.2 to the future of x — xi, is an apparent horizon, i.e., 
outgoing null rays have zero expansion on this hypersurface. Since a; = a;2 is also a singular point of the differential 
equations (2.11)-( ^.15 ), there is some concern that the hypersurface may be geometrically singular too. Nevertheless, 
since r is finite along x = X2 (except at u = 0) and g/g = 1 — 2m/r = 0, we conclude that the mass function and all 
curvature invariants are finite here, as are the scalar field and its derivatives (in a regular coordinate system). Thus, 
the apparent horizon is non-singular, and the region to its future is trapped. We have not continued the solution to 
larger values of x, but from the singularity theorems we know that to the future of the apparent horizon there is a 
central r = strong curvature singularity. This solution belongs in Class 11(a) described in Table I of Ref. [Tl| ]. 



D. Stability analysis 



Having constructed this CSS scalar field solution, we now turn to an examination of its stability properties. As 
mentioned in the introduction, we expect a critical solution to have a single unstable mode, an insight that results 
from efforts to understand the black hole mass scaling near the critical point 
the fields about its self-similar background value: 



13 1 . We therefore expand each of 
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■ 9o{x) +gi{x,T) , 
■9o{x)+gi{x,T) , 

-■jo{x) +jl{x,T) , 

: ho{x) + KT + hi{x, t) 



(2.22) 
(2.23) 
(2.24) 
(2.25) 



The subscript indicates a regular solution of Eqs. (2.11)-( p.l5 ); the subscript 1 is used to indicate linear perturbations 
about this solution. In terms of these fields, the perturbation equations become: 



2a; (ji 



X hi 

xg'i 

X9'i 
hiY 



■ Ji , 

91 - 9l : 

^ 47r [guo + 2.gojoii] , 
xiVo- 2x)j[ +ji{go 
+ jo9i +gixj'o ■ 



2x) 



(2.26) 
(2.27) 
(2.28) 

(2.29) 

To satisfy elementary fiatness at the origin, we set 51(0, r) = 5^(0, t) = and ji{0, t) — 0. 

We have used two independent numerical methods to search for the growing modes of Eqs. (2.26)-(2.29), a shooting 
method and a matrix eigenvalue method. Each method has particular strengths. Shooting allows accurate determina- 
tion of the unstable modes but relies on an initial guess close to the ultimate answer. For this reason, it is difficult to 
confirm that all the relevant modes have been found without further analysis. The matrix eigenvalue method appears 
to be less accurate than shooting, but it provides a scheme to determine all the modes at one time and thus confirm 
that we have identified all growing modes. The observed agreement between the two methods confirms that there 
is one unstable mode and two gauge modes. We describe each method in turn; Table || summarizes our numerical 
results. 



1. Shooting method 



This method provides direct computation of the cigcnmodcs of the perturbation equations (2.26)-(2.29) by recasting 
them as an eigenvalue problem. Consider solutions which separate in the form /i(a;,r) — e^'^ fi{x;X), where A is a 
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complex number, and the fields A) are also complex. Equations ( ^.26| )-(2.29) then reduce to a set of ordinary 
differential equations 



xh^= ji , 
xg'i — An 



25ojoji 



22;A(ji + hi) 



x{go - 2x)j[ + ji{ga - 2x) 
+ jogi+gixj'Q . 



(2.30) 
(2.31) 
(2.32) 

(2.33) 



Growing modes have real, positive A demonstrating the existence of an instability. As explained above, analyticity of 
the solutions at the origin (x = 0) and the similarity horizon [x — 1) is required since the critical solution observed 
in numerical simulations must be smooth everywhere. Analyticity determines the power series solution about a; = 
up to three real parameters corresponding to tan~^[(Imji)/(Reji)]|x=o, ReA and ImA. The amplitude |ji(0)| sets 
the scale of the perturbations and may be set to unity without loss of generality. Similarly, the power series solution 

about X = 1 is determined up to two complex (four real) parameters-ji(l) and hi{l). Equations (2.30)-( 2.33| ) were 
solved using a 4th order adaptive Runge-Kutta scheme, shooting from x = and x = 1 to a fitting point at x = 0.5 
until the seven parameters were adjusted to produce a smooth solution everywhere on the interval. 

There is one physical unstable mode with A ~ 1.065 when Attk^ = 0.577. This suggests that the critical exponent is 
7 = 1/A ~ 0.939, which is consistent with the numerical results for stiff fluid collapse where 7 < 0.96 In addition 
to the unstable mode, there are modes with A ~ 1.000 and A ~ 0.000. These modes correspond to the symmetries 
u u + const, and ip ip + const, of the background CSS solutions respectively, and do not change the physical 
spacetime. Thus, these are "gauge modes", and finding them at the right place is a (weak) test of our numerical 
methods. 

For completeness, we also searched for unstable modes of the CSS solution with Attk'^ = 1/3, but found only the 
two gauge modes as expected. 



2. Matrix eigenvalue method 



As the CSS solution presented above is not the threshold solution in massless scalar field critical phenomena, it 
is important to ensure that we have not missed any other unstable modes. One method to do this is to write the 
perturbation equations as u^r = Lu and to look directly for eigenvalues A of the time evolution operator L. 

For this, we have used the background and perturbation equations in polar-radial coordinates (they are given for 
example in [l4|] ), but this is accidental: we could equally well have implemented the time evolution operator method 
in Bondi coordinates, or the shooting method in polar-radial coordinates. 

The equations for the linear matter perturbation can formally be written as the system 

u^r = A{x)u^^ ^ B{x)u^ C{x)w, (2.34) 

where w(x,t) stands for the two first-order matter perturbation variables and w{x,t) stands for the two metric 
perturbations. (This is true both for Bondi coordinates and polar-radial coordinates.) The metric perturbations are 
not evolved, but are reconstructed from the matter perturbations by ODEs of the form 

w.-^=D{x)w + E{x)u. (2.35) 

The matrices A, B, C, D and E depend on the background solution. 

The perturbation equations were discretized in space but not in time. The resulting set of equations can formally be 
written as dujsf /dr = Lnun. Here un is the approximation of u on a grid with N points in the range < x < 1, and 
the matrix Ln is the corresponding finite difference approximation of the time evolution operator. The results quoted 
below were obtained at a resolution of A'' = 800, making a 1600^ matrix. The eigenvectors and eigenvalues of this 
matrix were then computed to machine precision using a standard linear algebra package. Most of the eigenvectors 
of Ln depend on the finite differencing scheme, but those with the highest values of ReA converge to modes of the 
continuum equations with increasing N, and the corresponding eigenvalues converge to the continuum eigenvalues A. 

The matrix eigenvalue method, in polar-radial coordinates, finds the growing mode at A ~ 1.055, and the gauge 
modes at A ~ 1.013 and A = 0. This last value is exact to machine precision as the ip ^ tp + const, gauge mode 
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decouples in these coordinates. Using the second-order convergence of our finite difference scheme, we estimate the 
uncertainty in these resuhs to be 0.01 for our finest grid of TV = 800 points. This estimate is consistent with the error 
in A for the gauge mode, which should be A = 1, as well as with the independent results from the shooting method 
in Bondi coordinates. 

With ijv already obtained, it was straightforward to implement the Lyapunov analysis described by Koike, Hara 
and Adachi and compare it with our matrix eigenvalue method. We just had to combine ijv with a fourth-order 
Runge-Kutta discretization in r in order to obtain a finite difference version of the evolution equations. (Such a 
method is variously called a method of lines, or a semi-discrete method.) We then used the Lyapunov method to 
obtain the highest few eigenvalues. In the limit in which the timestep Ar is chosen much smaller than the spatial 
resolution Ax, the resulting values of A coincide with those of the direct eigenvector method to machine precision, as 
one would expect. However, our method has the crucial advantage of producing all the top eigenvectors (modes) as 
well as the eigenvalues, whereas the Lyapunov method produces only the one eigenvector associated with the highest 
eigenvalue. Finding subdominant eigenvectors is important here because there are two gauge modes in the system of 
equations whose eigenvalues may be higher than that of the unstable physical mode. 



E. Comparison with the observed stiff fluid critical solution 



Having constructed a CSS scalar field solution with a single unstable mode, we now compare it with the stiff fluid 
critical solution. To facilitate the comparison, we adopt the polar-radial coordinates, {t,r), and variables (fluid and 
metric) used in the fluid studies H,^. This is effected by introducing t = — r^(a;) where 

x{\nO' = ^ , (2.36) 
g-x 

and ^(1) = 1. The metric and fluid variables are 

a' = g/g , (2.37) 

Lu = 4:nr^a^p^2TT[2x{f + K)/g- f] , (2.38) 

Fig. H shows our CSS scalar field solution and the stiff perfect fluid critical solution computed by evolving tuned 
initial data, as described in Ref. jl^,^. There is excellent agreement between the scalar field solution and the fluid 
solution inside the similarity horizon (^ = 1), and slightly beyond. The two solutions differ when ^ > ~ 1.14, 
however. Here the energy density, p, of a fluid equivalent to the scalar field vanishes and the fiuid 3-velocity reaches 
unity. In the numerical evolutions of the perfect fiuid, a floor is imposed to maintain p > and a timelike 4- velocity in 
the simulations. Th is acc ounts for the dramatic difference between the solutions at larger values of ^, and is discussed 



in detail in Section [II B 



In addition to this direct comparison of the perfect fluid and scalar field solutions, we verify that the solutions share 
the same stability properties. As mentioned previously, the unstable mode of the CSS scalar field solution gives a 
critical scaling exponent of 7 ~ 0.94, which agrees with the perfect fiuid estimate of 7 < 0.96. Finally, we can directly 
compare the unstable modes as shown in Fig. 0. 



III. THE ROLE OF THE SCALAR FIELD CSS SOLUTION IN GRAVITATIONAL COLLAPSE 



There are at least three regular self-similar, spherically symmetric scalar field solutions: the cosmological CSS 
solution found by Christodoulou, the CSS solution with one growing mode found in stiff ffuid collapse and studied 
here, and the DSS critical solution with one growing mode ^ examined in more detail in Ref. Turning to the 
questions raised in the Introduction, we examine why the stiff ffuid and scalar field critical solutions are not identical, 
given the (local) equivalence between the two matter models. 



A. Self-similar solutions in massless scalar field collapse 



The cosmological solution [Attk^ — 1/3) has additional symmetries; it is actually a spatially fiat Friedmann solution 
driven by an homogeneous scalar field. It has no growing perturbations and is therefore a global attractor. The 
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FIG. 2. The CSS solution with 4ttk.^ ~ 0.577 (solid lines), obtained by solving the ODEs in Eqs. ( ^.17| )-( |2!2l| ) and Eq. ( ^.36| ), 
compared to a near critical solution computed using the PDE code in Ref. (dotted lines). The similarity horizon is at ^ = 1. 
There is excellent agreement between the two solutions up to, and sightly beyond, the similarity horizon. The self-similar scalar 
field solution cannot be identified with a stifT fluid solution beyond ^ ~ 1.14 where the gradient of the field is spacelike. 



singularity in this cosmological solution is velocity dominated. Thus, as one approaches the spacelike singularity, the 
evolution of each point in space decouples from that of neighboring points and evolves approximately as a Friedmann 
solution In general, we expect this solution to be the final attractor inside the black hole as the central singularity 
is approached. 

The scalar field CSS solution with Attk^ ~ 0.577 has one growing mode, and therefore, one might expect it to be a 
scalar field critical solution. However, a critical solution is defined through its appearance at a black hole threshold, 
and thus the solution's global structure must also be considered. Examining the solution beyond xi, where the scalar 
field gradient becomes spacelike, we find that this solution has an apparent horizon at some x = X2, where the mass, 
TO, and radius, r, satisfy to = r/2. This means that the candidate critical solution is not actually on the black hole 
threshold, and that is why it is not a critical solution of scalar field collapse even though it has precisely one growing 
mode. 



B. Self-similar solutions in stiff fluid collapse 



We now consider the scalar field CSS solutions discussed above as candidate solutions for stiff fluid critical collapse. 
We can quickly dispatch with the DSS scalar field solution as a possibility. As shown in Fig. ^, the equivalent stiff 
fluid density, p, in the DSS solution does not have a definite sign; it may be positive, zero, or negative for some r at 
any x. Therefore, the scalar field cannot be interpreted globally as a fluid. As mentioned in the introduction, this has 
been well known and understood for some time. 

In the previous section, we found that the Attk^ ~ 0.577 CSS solution matches the intermediate attractor in stiff 
fluid collapse, both in functional form and in the stability properties. However, this solution is not at the threshold 
of black hole formation in scalar field collapse, which raises the immediate question of why it is at the black hole 
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FIG. 3. This figure sliows tlie unstable mode for the metric function a, comparing the mode of the CSS scalar field solution 
(solid line) with that of the perfect fluid (triangles). The modes match particularly well when ^ < (^i — 1.14), where the 
two solutions are equivalent. Beyond this point, where the fluid and scalar fleld are no longer equivalent, the modes diverge as 
expected. The perfect fluid mode was calculated using from two near-critical fluid simulations. As the mode has no intrinsic 
magnitude, the fluid data were scaled to match the scalar field mode. 




FIG. 4. The equivalent fluid density, given by Eq. (2.3), of the DSS scalar field solution, plotted for < r < A ~ 3.44 at 
an arbitrarily chosen point, in this case x = 0.5. 



threshold as a stiff fluid. The answer to this question comes in demanding that the fluid solution (putatively equivalent 
to the scalar field) be physical. 

Consider the CSS perfect fluid critical solutions obtained from solving a CSS ansatz for different values of k. (See 
Hara, Koike and Adachi Q for a derivation of the equations, and a solution method. See |^ for the solutions with 
k > 0.89.) Figure ^ shows the fluid density uj{£^) as k ^ 1. For fc < 1, w is smooth, monotonically decreasing, and 
asymptotes to zero from above. As /c — s- 1 from below, the curve develops a corner joining two approximately 
linear sections. In the limit k = 1 one would expect w(^) to reach zero at finite ^ with finite nonzero slope, and to 
continue as w = for ^ > ^i, with the corner becoming a kink at ^i. But this is not the case. For £,<£,!, the scalar 
field solution is indeed the limit fc ^ 1 of the fc < 1 fiuid solutions, but for ^ > it is qualitatively different: it 
continues smoothly through to negative values of uj. 

Mathematically, it is a matter of definition how one continues stiff fluid solutions to the future of points of zero 
density. One choice is to evolve the equivalent scalar field, which remains well-defined through this transition region. 
The density defined in this way becomes negative but is smooth everywhere. The 4- velocity is not smooth where the 
density changes sign, but the three-velocity v with respect to constant r observers is: it continues smoothly to values 
|u| > 1. A second choice, motivated by the desire for a "physical" fluid solution, is to truncate the CSS solution 
at X = xi, and replace the spacetime outside this point with another solution of the Einstein equations. It is not 
possible, however, to match to a vacuum Schwarzschild solution without invoking some sort of thin shell since the T^^ 
component of the stress-energy tensor does not vanish at xi although the (equivalent) fluid energy density does. The 
divergence of the Lorentz factor, (1 — w^)~^/^, at xi keeps finite as p — > 0. 

In the numerical simulations ^ of stiff fluid critical collapse, the fluid solution is also modified for a; > xi. During 
the evolution the code enforces the physical fiuid conditions by imposing a fioor, or minimum relative value of the 
fiuid energy density to the momentum, on the solution. (See Appendix ^ for information regarding the fioor, and 
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FIG. 5. A detail of the region about showing uj for the perfect fluid CSS solutions (obtained from a CSS ansatz), including 
fc < 1. 

Ref. jl^ for additional information on this code.) The floor is a common, ad hoc method for allowing vacuum or 
near-vacuum regions to develop in numerical solutions. Examining a near-critical, stiff fluid evolution, we flnd that 
the floor is only active for x>^ xi. One visible effect of the floor on this solution is the slight difference in the evolved 
fluid solution with the ODE solution near a; < cci, as shown in Fig. ^j. This modiflcation, however, does not change 
the scaling properties of the critical solution. The eigenvalue spectrum of a self-similar solution depends only on the 
solution inside the past light cone of the singularity, i.e., the region x < xi for the stiff fluid solution where the floor 
is not active. 

However the CSS solution for the stiff ffuid is modified beyond xi, the effect is the disappearance of the apparent 
horizon observed in the scalar field solution. Thus, this solution can appear at the threshold of black hole formation 
in stiff ffuid spacetimes, and is indeed the observed critical solution. This view is further bolstered by the apparent 
continuity of the critical solutions in /c-space for fluid collapse with equations of state of the form P = kp, 0.05 < k < 
1 0- 

IV. CONCLUSIONS AND DISCUSSION 

Scalar fields and perfect fluids have historically been important models in establishing our understanding of critical 
phenomena in gravitational collapse. The dominant features of critical behavior are now understood for both models 
in spherical symmetry. The real scalar field critical solution is discretely self-similar. Type II, and has a mass-scaling 
exponent 7 ~ 0.37. The perfect fluid critical solutions are also Type II, continuously self-similar, and the stiff (P = p) 
fluid's mass-scaling exponent measured in near-critical simulations is 7 < 0.96. These models also share a well known 
formal equivalence, at least locally, between irrotational stiff fluids and real scalar flelds. Thus, the very different 
critical behavior observed in these two possibly equivalent models poses an interesting question: Why is a critical 
solution for one model not the critical solution in the other? 

At least part of the answer comes in how one defines a fluid. In the usual construction of the fluid model as 
the continuum limit for an underlying discrete system, the fluid satisfles certain physical conditions, e.g., that the 
co-moving energy density is positive, and the four-velocity timelike. These conditions limit the scalar fleld solutions 
that can be interpreted as fluids. One such disqualified solution is indeed the DSS scalar field critical solution, for 
which the equivalent fiuid density becomes negative. 

To better understand the stiff ffuid critical solution, we have explicitly constructed it as a boundary value problem, 
using the scalar field variables to describe the fiuid. This confirms that an exactly CSS, regular, stiff-ffuid solution 
exists. It has exactly one growing perturbation mode, whose inverse Lyapunov exponent 7 ~ 0.94 is equal, to within 
estimated numerical error, to the critical exponent computed from critical collapse of a stiff fluid. This CSS solution 
is not seen as a critical solution in scalar field collapse because it contains an apparent horizon, and is therefore not 
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at the black-hole threshold. 

In spite of the fact that the complete scalar field solution contains an apparent horizon, and it cannot be interpreted 
globally as a fluid, this CSS solution does appear in critical collapse of the stiff fluid. The apparent horizon of the 
CSS solution is in a region of the critical solution where the scalar field gradient is spacelike, so that the equivalent 
fluid density is negative. However, the solution allows a physical fluid interpretation inside the similarity horizon, and 
the collapsing fluid apparently is attracted to this part of the CSS solution. 

Generic stiff-fluid solutions reach points of zero density and lightlike 4-velocity in a finite time, at which point the 
fluid equations break down. One can either continue the solution as a massless scalar field or one can replace the 
unphysical fluid with another solution of the Einstein equations. In computational simulations of fluid critical collapse, 
the solutions are modified beyond x = xi via a floor to maintain a positive fluid density. The resulting near-critical 
spacetimes form a continuous family parameterized by k. The k — 1 solution approximates the CSS solution until 
just before the spacelike surface (x — x\ in our coordinates) where the density goes to zero. Beyond that surface, 
they are qualitatively different — matter at a; > a:i is ejected at almost the speed of light, and no apparent horizon 
forms. This allows the modified CSS solution to function as a critical solution at the black hole threshold. 

Finally, we note the recent work of Harada ||2C|l , who argues that an additional unstable mode (the so called 
"kink mode"), is present in self-similar fluid solutions with P = kp, for k > 0.89. This mode is characterised by a 
discontinuity in the derivative of the fluid density (taken with respect to the similarity variable) at the sonic point. 
To date, there has been no evidence in near-critical fluid evolutions of such additional unstable modes, although if the 
growth factors are sufflciently small, it is entirely possible that they would simply not be visible given the numerical 
precision typically used in the PDE evolutions. In the context of the current work, our requirement of analyticity at 
the sonic point eliminates any possibility of our finding a kink mode via the CSS ansatz. Further work will be needed 
to clarify this situation. 
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APPENDIX A: EQUIVALENCE BETWEEN AN ARBITRARY FLUID AND A SCALAR FIELD 

An irrotational perfect fluid is equivalent to a scalar field for any one-parameter equation of state P — Pip) |ll8| , |l9| . 
The conservation of the perfect fiuid stress-energy tensor 

Tab = (-P + p)UaUb + Pgab (Al) 

is equivalent to the 3-f 1 equations 

(P + p)u"Va1i^ + h'^'^WaP = 0, (A2) 

u-^VaP + (P + p)VaM° = 0. (A3) 

Here hab = gab + UaUb is the projector into the 3-space orthogonal to the fiuid 4-velocity, it". The first equation is a 
3-vector equation, and is commonly known as the force, or Euler equation. The second equation is a scalar equation 
expressing the conservation of mass/energy. We now parameterize p and it" together through a single vector that is 
not a unit vector: 

w" = h{p)u'' ^ h{pf = -Waw"-. (A4) 
If we now choose the function h{p) to be a solution of the ordinary differential equation 

dJ^_JP_ 

h - P + p' ^^^> 
then we find that the Euler equation takes the simple form 

w'^h^'^VaWc - VcWa) = 0. (A6) 
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(Up to a factor, ( [A5| ) defines h to be the enthalpy per particle of the fluid.) If the fluid is irrotational 

V[aU6] = 0, (A7) 



we also have h'^'^h^'^V^a'Wi,] = 0. Together with (A6), we then have Vja^b] = 0, and therefore, Wa can locally be 
written as the gradient of a scalar function, Wa — Va0- Expressed in terms of alone, the remaining field equation 
(A3) becomes 

Here }\} = —W^ipWacj), and dp/dP is a given function of ft,, determined by the equation of state P = P{p) and the 
solution h — h{p) of ([A5|). Using m° = h^^V^cj) and defining the projector /i*^^ as before, and with the sound speed 
squared given by = dP/dp, the characteristics of the scalar field equation can be emphasized by writing it as 

(-uV + c^/i'^^) VaVh(/) = 0. (A9) 



APPENDIX B: STIFF FLUID AND SCALAR FIELD VARIABLES 



Fluids are usually described analytically using the fundamental quantities of density, p and four- velocity, Modern 
numerical methods for solving the fluid equations, however, are designed to exploit conservation properties, and thus 
describe the fluid in terms of its conserved properties, e.g., the fluid's energy and momentum. Modelling a dynamic 
stiff fluid can be challenging, and a new set of conservation variables 

Ilfiuid = {au^'fp + {au^-fp + (au'~)(au*)(p + p), 

$fluid = {avJ'fp + {av^ fp - {au'){av'){p + p). (Bl) 

was found that considerably improves numerical stability and performance jlTt . These variables are linear combina- 
tions of components in the stress-energy tensor, representing energy and momentum, and both are positive-definite 
when the fluid 4-velocity is timelike, and the density and pressure are positive. Nevertheless, numerical error in a 
computer evolution will sometimes result in values for <&fluid and/or Ilfluid that are negative, which translates into a 
negative fluid density. In order to preserve the correct physics, a "floor" is imposed on these variables by forcing them 
to always be greater than a small, positive constant, 5. This floor is applied at every step in the fluid update as 

Ilfluid <^ max((5, Ilfluid), (B2) 
^fluid ^ max((5, <I>fluid)- (B3) 

5 must be small enough that its influence on the dynamics of the physical fluid is minimized, and its effect can, at 
least partially, be judged by repeating the simulations of the same initial data while varying 5. In the critical collapse 
study of Ref. (|], the floor was set to 5 = 10~^°. 

Expressing these new fluid conservation variables in terms of a set of first-order variables for the scalar field provides 
an interesting comparison. A natural choice of such variables [0 is 

^scalar — V,r, Hgcalar = —^.t- (B4) 

a 

These are related to the fluid variables for the stiff fluid {P — p) as 

Ilfluid = 2^ ^(^scalar + Hgcalar)^, 

^'fluid = ^a"^($scalar - Uscalar)^- (B5) 

We see that the fluid variables are essentially squares of the scalar field variables. The linear scalar field equation 
(written in first order form) is equivalent to the nonlinear fiuid equations, and can be recovered by changing to the 
"square-root variables" . The fluid variables are positive by definition, but one of them goes to zero when the scalar 
field variables change sign, as 

(0^^0'^)2 = fflfluid^fluid (B6) 

Keeping the fluid variables away from zero prevents such a sign change and qualitatively changes the solutions. 
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